home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YDES.C < prev    next >
Text File  |  1988-11-28  |  7KB  |  330 lines
  1. /******************* (C) 1986,7,8 by JAMES A. YORKE **************************/
  2. /******************************* YDEs.C ***************************************
  3.     This file is where the differential equations are defined.  They are
  4.     called (generally) by DifEqu() which is called by rungekutta(),
  5.     and rungekutta() is in YTOOLS.C.   Notice that Y[] is a temporary 
  6.     position vector which is set by rungekutta to be either y[] or temp[].
  7. DED(k,Y)  DUFFING: x'' + C1*x' -C2*x +C3*x*x*x = rho*sin(C4*t) 
  8. initDED()
  9. DEPD(k,Y) PARAMETRICALLY DRIVEN DUFFING: x'' + C1*x' 
  10.    -x + C2*(1+C5*sin(C4*t))*x +C3*(1+C6*sin(C4*t))*x*x*x = 0 
  11. initDED()
  12. DEP(k,Y) Forced Damped Pendulum x''+C1*x'+C2*sinx =rho*(C3+cos(C4*t))
  13. initDEP()
  14. DEV(k,Y)  V:  a VanDerPol differential equation due to UEDA 
  15. initDEV() for V:  a VanDerPol differential equation due to UEDA 
  16.  
  17. ********************** (C) 1986,7,8 by JAMES A. YORKE ************************/
  18. #include "yinclud.h"
  19.  
  20.  
  21.  
  22.  
  23.  
  24. /*  Y varies in different stages of the rungekutta 
  25.         and is usually not equal to y */
  26.  
  27. DED(k,Y) /* DUFFING: x'' + C1*x' -C2*x +C3*x*x*x = rho*sin(C4*t) */
  28.     double k[],Y[];
  29. {
  30.     int base,i;
  31.     double t,x,y,xx,xxx;
  32.     t = Y[0];
  33.     x = Y[1];
  34.     y = Y[2];
  35.     xx = x*x;
  36.     xxx = x*xx;
  37.  
  38.     k[0] = 1;
  39.     k[1] = y ;
  40.     k[2] = -C1*y + C2*x -C3*xxx + rho*sin(C4*t) ;
  41.     k[3] = 1;  /* coordinate 0 is set to 0 after each cycle */
  42.      if (num_lyap > 0) 
  43.      {
  44.     for (i = 0; i < num_lyap; i++)
  45.     { 
  46.       base =  lyapzero + vec_dim*i;
  47.       k[base] =     Y[ base + 1];
  48.       k[base +1 ]  = (C2 - 3*C3*xx) *Y[ base ] -C1*Y[ base + 1];
  49.     }
  50.      }
  51. }
  52. initDED()
  53. {
  54.     int d_mod();
  55.     zeroth = 1;
  56.     vec_dim = 2;   /*the dimension of the 
  57.         Lyapunov vectors = phase space dim*/
  58.     num_lyap = 0; /*the number of Lyapunov numbers 
  59.         to be computed <= vec_dim */
  60.     lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
  61.         zeroth lyapunov vector */
  62.         /* variables: time mod twopi/C4 ,x,y,t */
  63.  
  64.     X_upper = 2.2 ;    /* x scale */
  65.     X_lower = -2.2;
  66.     Y_upper = 2.;     /* y scale */
  67.     Y_lower = -2.;
  68.     X_coord = 1;
  69.     Y_coord = 2;
  70.     C1 = 1.;
  71.     C2 = 10;
  72.     C3 = 100;
  73.     rho = .85;
  74.     C4 = 3.5;
  75.     modPointer = d_mod;
  76.     DE = DED;   /*  DE is a pointer to a function   */
  77.     steps_per_cycle = 80;
  78.     /* now PickMap in YPickMap.C will automatically set:
  79.             step = (twopi/C4)/steps_per_cycle; */
  80. }
  81. d_mod()/* this is invoked by Duffing and by Parametrically Driven Duffing */
  82. {
  83.     double moduloAB();
  84.     modFlag = NO;
  85.     y[0] = moduloAB(y[0],0.,twopi/C4);
  86. }
  87.  
  88. /* DEPD PARAMETRICALLY DRIVEN DUFFING: 
  89.    x'' +C1*x' -x +C2*(1+rho*sin(C4*t))*x^2 +C3*(1+rho*sin(C4*t))*x^3 = 0 */
  90. DEPD(k,Y)
  91.     double k[],Y[];
  92. {
  93.     int base,i;
  94.     double t,x,y,xx,xxx,s,xxCoef,xxxCoef;
  95.     t = Y[0];
  96.     x = Y[1];
  97.     y = Y[2];
  98.     xx = x*x;
  99.     xxx = x*xx;
  100.     s = 1 + rho*sin(C4*t);
  101.     xxCoef = C2*s; 
  102.     xxxCoef = C3*s;
  103.     k[0] = 1;
  104.     k[1] = y ;
  105.     k[2] = -C1*y +x -xxCoef*xx -xxxCoef*xxx;
  106.     k[3] = 1;  /* coordinate 0 is set to 0 after each cycle */
  107.      if (num_lyap > 0) 
  108.      {
  109.     for (i = 0; i < num_lyap; i++)
  110.     { 
  111.       base =  lyapzero + vec_dim*i;
  112.       k[base] =     Y[base + 1];
  113.       k[base +1 ]  = (1. -2.*xxCoef*x - 3.*xxxCoef*xx) *Y[ base ] 
  114.         -C1*Y[ base + 1];
  115.     }
  116.      }
  117. }
  118.  
  119. initDEPD()/* Parametrically Driven Duffing */
  120. {
  121.     int d_mod();
  122.     zeroth = 1;
  123.     vec_dim = 2;   /*the dimension of the 
  124.         Lyapunov vectors = phase space dim*/
  125.     num_lyap = 0; /*the number of Lyapunov numbers 
  126.         to be computed <= vec_dim */
  127.     lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
  128.         zeroth lyapunov vector */
  129.         /* variables: time mod twopi/C4 ,x,y,t */
  130.  
  131.     X_lower = -5.;
  132.     X_upper = 3. ;    /* x scale */
  133.     Y_lower = -3.;
  134.     Y_upper = 6.;     /* y scale */
  135.     X_coord = 1;
  136.     Y_coord = 2;
  137.     C1 = 1.;
  138.     C2 = 1.;
  139.     C3 = 1.;
  140.     C4 = 1.;
  141.     rho = .5;
  142.     rho_step = .05;
  143.  
  144.  
  145.     modPointer = d_mod;
  146.     DE = DEPD;   /*  DE is a pointer to a function   */
  147.     steps_per_cycle = 50;
  148.     /* now PickMap in YPickMap.C will automatically set:
  149.             step = (twopi/C4)/steps_per_cycle; */
  150. }
  151.  
  152. /*    y[0] = moduloAB(y[0],0.,1.);
  153.     y[1] = moduloAB(y[1],0.,1.);
  154. */
  155.  
  156.  
  157.  
  158.  
  159. DEP(k,Y)/* Forced Damped Pendulum x''+C1*x'+C2*sinx =rho*(C3+cos(C4*t))*/
  160.     double k[],Y[];
  161. {
  162.     int i;
  163.     double t,x,y,sinx,C2cosx;
  164.     int base;
  165.     t = Y[0];
  166.     x = Y[1];
  167.     y = Y[2];
  168.     sinx=sin(x);
  169.  
  170.     k[0] = 1;
  171.     k[1] = y;
  172.     k[2] = -C1*y -C2*sinx + rho*(C3+cos(C4*t));
  173.     k[3] = 1;  /* coordinate 0 is set to 0 after each cycle */
  174.  
  175.      if (num_lyap > 0) 
  176.      {
  177.     C2cosx = C2*cos(x);
  178.  
  179.     for (i = 0; i < num_lyap; i++)
  180.     { 
  181.       base =  lyapzero + vec_dim*i;
  182.       k[base]
  183.        =     Y[ base + 1];
  184.       k[base +1 ] 
  185.        = -C2cosx *Y[ base ]
  186.              -C1*Y[ base + 1];
  187.     }
  188.      }
  189. }
  190.  
  191.  
  192. /*    y[0] = moduloAB(y[0],0.,1.);
  193.     y[1] = moduloAB(y[1],0.,1.);
  194. */
  195.  
  196. v_mod()
  197. {
  198.     double moduloAB();
  199.     y[0] = moduloAB(y[0],0.,twopi/C4);
  200.     modFlag = NO;
  201. }
  202. p_mod()
  203. {
  204.     double moduloAB();
  205.     y[0] = moduloAB(y[0],0.,twopi/C4);/* time */
  206.     modFlag = NO;
  207.     if (X_coord != 1)
  208.         y[1] = moduloAB(y[1],-pi,pi);
  209.     else 
  210.         y[1] = moduloAB(y[1],X_lower,X_lower + twopi);
  211. }
  212.  
  213.  
  214. y_mod_for_des()
  215. {
  216.     return;
  217. }
  218.  
  219.  
  220. initDEP()
  221. {
  222.     int mod_p();
  223.     zeroth = 1;
  224.     vec_dim = 2;   /*the dimension of the 
  225.         Lyapunov vectors = phase space dim*/
  226.     num_lyap = 0; /*the number of Lyapunov numbers 
  227.         to be computed <= vec_dim */
  228.     lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
  229.         zeroth lyapunov vector */
  230.         /* variables: time mod twopi/C4 ,x,y,t */
  231.  
  232.     X_upper = pi ;    /* x scale */
  233.     X_lower = -pi;
  234.     Y_upper = 4.;     /* y scale */
  235.     Y_lower = -2.;
  236.     X_coord = 1;
  237.     Y_coord = 2;
  238.     C1 = .2;
  239.     C2 = 1.0;
  240.     C3 = 0.;
  241.     rho = 2.5;
  242.     C4 = 1.;
  243.     steps_per_cycle = 30;
  244.     period = 1;/* for Newton method */
  245.     its_per_plot = steps_per_cycle;
  246.     /* now PickMap in YPickMap.C will automatically set:
  247.             step = (twopi/C4)/steps_per_cycle; */
  248.     DE = DEP;   /*  DE is a pointer to a function   */
  249.     modPointer = p_mod;/* takes angle mod two pi, more or less */
  250. }
  251.  
  252.  
  253. DEV(k,Y) /* V:  a VanDerPol differential equation due to UEDA */
  254.     double k[],Y[];
  255. {
  256.     int i,base;
  257.     double t,xx,x,y;
  258.     t = Y[0];
  259.     x = Y[1];
  260.     y = Y[2];
  261.     xx = x*x;
  262.     k[0] = 1;
  263.     k[1] = y;
  264.     k[2] = C1*y*(1-xx) -C2*x -C3*x*xx + rho*sin(C4*t);
  265.     k[3] = 1;  /* coordinate 0 is set to 0 after each cycle */
  266.      if (num_lyap > 0) 
  267.      {
  268.     for (i = 0; i < num_lyap; i++)
  269.     { 
  270.       base =  lyapzero + vec_dim*i;
  271.       k[base] =     Y[ base + 1];
  272.       k[base +1 ]  = (-2*x*C1*y -C2 - 3*C3*xx) *Y[ base ] 
  273.             +C1*(1-xx)*Y[ base + 1];
  274.     }
  275.      }
  276. }
  277.  
  278. initDEV()/* for V:  a VanDerPol differential equation due to UEDA */
  279. {
  280.     int v_mod();
  281.  
  282.     zeroth = 1;
  283.     vec_dim = 2;   /*the dimension of the 
  284.         Lyapunov vectors = phase space dim*/
  285.     num_lyap = 0; /*the number of Lyapunov numbers 
  286.         to be computed <= vec_dim */
  287.     lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
  288.         zeroth lyapunov vector */
  289.         /* variables: time mod twopi/C4 ,x,y,t */
  290.     dim = lyapzero + num_lyap*vec_dim;
  291.             /* needed for rungekutta */
  292.  
  293.     X_upper = 2. ;    /* x scale */
  294.     X_lower = -2.;
  295.     Y_upper = 2.;     /* y scale */
  296.     Y_lower = -2.;
  297.     X_coord = 1;
  298.     Y_coord = 2;
  299.     C1 = .2;
  300.     C2 = 0.;
  301.     C3 = 1.;
  302.     rho = 2.0;
  303.     C4 = 1.;
  304.     rho = 2.0;
  305.     steps_per_cycle = 40;
  306.     /* now PickMap in YPickMap.C will automatically set:
  307.             step = (twopi/C4)/steps_per_cycle; */
  308.     DE = DEV;   /*  DE is a pointer to a function   */
  309.     modPointer = v_mod;/* specifies a modulo calculation */
  310. }
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.